Surface premelting and melting of colloidal glasses

The nature of liquid-to-glass transition is a major puzzle in science. A similar challenge exists in glass-to-liquid transition, i.e., glass melting, especially for the poorly investigated surface effects. Here, we assemble colloidal glasses by vapor deposition and melt them by tuning particle attractions. The structural and dynamic parameters saturate at different depths, which define a surface liquid layer and an intermediate glassy layer. The power-law growth of both layers and melting front behaviors at different heating rates are similar to crystal premelting and melting, suggesting that premelting and melting can be generalized to amorphous solids. The measured single-particle kinetics reveal various features and confirm theoretical predictions for glass surface layer.


Glass stability
Ultrastable glasses have been fabricated in metallic (80), molecular (32) and polymer (81) systems by vapor deposition, but not yet in colloids. Here, by increasing the sample temperature, i.e. enhancing the attraction, the colloidal liquid region starts to vitrify into a glass and continues to grow through the slow layer-by-layer vapor deposition for more than 10 h until the vapor phase is nearly depleted ( fig. S2 and movie S1). Such glasses formed by vapor deposition are highly stable due to effective surface relaxation (32,82). χ d = v d τ α s /σ is inversely proportional to the degree of surface relaxation which has been proved to be the dominant factor on glass stability (82). v d is the deposition rate, τ α s is the relaxation time of the surface layer, and σ is the average particle diameter. A glass formation process with a low χ d produces a highly stable glass. The time t 0 = σ 2 /(2(d s − 1)D) is the mean time for a particle on the vapor-liquid interface moving one mean diameter, and is commonly used as τ α s (83). D is the diffusion coefficient of a particle on the surface (84) Besides the surface relaxation rate, the bulk density of our glasses 0.83 is higher than typical colloidal glasses and close to the upper limit 0.855 of the randomly close-packed hard disks (88).
Moreover, the bulk elastic moduli in fig. S14 are one order of magnitude higher than those of normal colloidal glasses (89,90) and those of colloidal crystals of the same system (23), see the detailed calculation in Section B.5. In summary, the small χ d , melting from surface instead of bulk, and high elastic moduli suggest that our glasses appear to be relatively stable. On the other hand, it is challenging to measure the ultrastability directly from the relaxation time τ due to our limited experimental time scale.

Image analysis for large and small particles
Particles can be identified by conventional image analysis (33). To further distinguish large and small spheres, we measure the radius of gyration of the pixel brightness R g and the total brightness of each particle's image. These two parameters can distinguish large and small particles ( fig. S5).

Dividing the field of view into stripes
The field of view is divided into stripes parallel to the glass surface with a width of 30 pixel ≈ 1 σ ( fig. S7). We measure various quantities for each particle i and average them in each stripe region to obtain their profiles along the y direction.

Voronoi tessellation
We use the radical Voronoi tessellation in the Voro++ library (91) instead of traditional Voronoi tessellation for our binary system, because otherwise, the bisecting line between large and small spheres may cut through the large sphere and cannot reflect the real Voronoi area associated with each sphere.
In the 2D radical Voronoi tessellation, the radical curve is composed of the points with the same tangent length for the two neighboring spheres, thus the tangential lines from the radical curve to each sphere's surface have the same length. Radical Voronoi tessellation avoids the intersection through spheres and possesses the topological features of the tessellation. Fig. S8 illustrates an example of the radical Voronoi tessellation of our binary system. The local density of particle i, where A i is the area of its Voronoi polygon.

Overlap function
Overlap function has been used to define the glass melting interface in simulation (8). It reflects the similarity of particles' configurations after a time interval t (8): where Θ(x) is the Heaviside step function, N y is the number of particles in the stripe at y. If a particle moves more than a = 1.5 σ, it is deemed mobile and no longer contributes to F m even if it eventually returns close to its original position (8). Fig. S12A shows F m (t) under the fast temperature change.
The melting interface is defined as F m = e −1 , as shown in fig. S12A (8). Interestingly, the melting interface at every instance determined from the dynamic parameter F m (t, y) (red cross in fig. 12B) coincides with the interface between surface liquid and glassy layer defined by the structural parameter ρ = 95% in the main text.

Melting and glass transition temperatures
The glass transition temperature featured with infinite relaxation time is often identified as the modecoupling critical temperature T MCT measured by extrapolating the relaxation time to infinity using Eq. 5 from mode coupling theory (56,92). The fitted ρ c = 0.80 in Fig. 5 corresponds to the bulk density at 25.4 which is close to the melting temperature 25.3 for the slow temperature change.
Glass is also empirically defined when the material has become too viscous to flow in a reasonable time scale (56), e.g. the viscosity (or relaxation time) > 10 15 times of that of a normal liquid (56).
Given that viscosity changes rapidly near the glass transition temperature, T g is not sensitive to the choice of threshold value. Besides viscosity, properties such as heat capacity, density and elastic moduli suddenly change during cooling or heating at a certain rate (2); either the onset, the middle or the ending of such a change has been used as an empirical definition of T g (17). In the main text, T g is fitted from Eq. 1 on the penetration depth, which avoids distinguishing supercooled liquid and glass.
The fitted T g agrees with the abrupt change in bulk density, the extrapolated zero elastic moduli shown in fig. S14 and the fully melted sample from the direction observation (movie S2).
The elastic moduli in the main text are calculated according to the following four steps (23,93): (1) Construct the covariance matrix of particle displacements and calculate its eigenvalues and where a particle's coordinate u i ∈ x 1 , y 1 , . . . , x N , y N . t is the average over frames. i, j = 1, 2, . . . , 2N . N is the number of particles. We track the positions of N ≈ 10 3 particles in the bulk region for 12000 frames at 20 frames/s for measuring the high-frequency modes. We diagonalize the large 2N × 2N matrix C ij to solve the eigenvectors, which are the polarization vectors of the normal modes. The eigenvalues λ = k B T /(mω 2 ), which provide the angular frequencies ω of the corresponding modes.
The mass of a small particle m is mass unit. The accurate value of the particle mass is not needed because it is canceled in the final result. The displacement of a particle is in the unit of the average diameter σ. The wave vector q is in the unit of 1/σ, and ω is in the unit of k B T /mσ 2 . (2) Extract the dispersion relation ω(q) (fig. S13B). The Fourier decomposition of the eigenmodes into transverse and longitudinal components yields two spectral functions, E T (q, ω) and E L (q, ω), respectively: where q = q/| q|, is average over different directions of q, e j (i) = (e ix , e iy ) ω j is the polarization vector of the ith particle corresponding to the jth mode. ω(q) is obtained by selecting an ω that maximizes S13A). As at low q, ω and q exhibit a linear relationship (fig. S13B), thus we can obtain the longitudinal modulus, M = ρ 2D (lim q→0 (∂ω L /∂q)) 2 , the shear modulus G = ρ 2D (lim q→0 (∂ω T /∂q)) 2 , and the bulk modulus, B = M-G. The areal density ρ 2D = 6mh/(πσ 3 ) and the height of the sample cell h ≈ σ. The bulk and shear moduli are (38000, 18000) k B T/σ 2 for the vapor-deposited glass, which are much higher than (450, 100) k B T/σ 2 for the colloidal glasses (89,90) and (800, 400) k B T/σ 2 for the colloidal crystals (23). These results are converted into the same unit for comparison and are in accordance with the fact that stable glasses have larger elastic moduli.

Multilayer glasses
For bilayer and trilayer glasses, the bright-field images are blurry because the different refractive indexes of water and PMMA spheres produce diffractions from other layers in the z direction. We find that the dense region is brighter and the pixel brightness fluctuates more at regions with more mobile particles ( fig. S17, A and B). Therefore, the pixel brightness (B) and its fluctuation (i.e. the standard deviation of B(t) over 30 s, std(B)) can be used to quantify the coarse-grained local density and dynamics, respectively. This result is confirmed by the calibration in the monolayer sample in fig. S18.
We measure the mean brightness B of each coarse-grained 30×30 pixel 2 region and its log(std(B)) and then rescale them to [0, 1], denoted asB andlog(std(B)), respectively. In monolayer samples, their profiles well agree with the profiles ofρ andlog(DW ) respectively (fig. S18, A and B) and produce the same fitted T g and power law exponent α from d 1,2 ( fig. S18C). Similarly, we measure the coarse-grained

Landau theory
For a 2D vapor-solid system with a 1D free surface along the x direction, the general form of the Landau free energy functional (45) is where y s is the position of the vapor interface, and ψ is a scaler order parameter such as density. The three terms in Eq. S4 represent the interface, the surface region with a gradient and the homogenous bulk region, respectively. The Taylor expansion of the bulk term f (ψ) is (45) f Minimising F {ψ} by δF/δψ = 0 yields (45) Eq. S6 has been solved analytically in ref. (45), and we rewrite the solution in the following form (45): where Q, P and S are parameters related to temperature.
is the bulk value at the melting point; is the bulk value. ψ b can be directly measured; therefore, the three free parameters a, b and ψ(y = y s ) can be solved from the fitted Q, P and S from Eq. S7.
Substituting Eq. S6 into Eq. S4, we obtain the free energy of the whole system in Eq. S4 (45): Because the last term ∞ ys f (ψ b )dy is a constant, we only need to consider the first two terms (45) Under a given vapor interface position y s , the fitting of Eq. S7 gives a, b, c in Eq. S5. By tuning y s , we obtain the corresponding F in fig. S21A. The minimum of F gives y s and ψ(y) of the equilibrium state. The bulk energy is shown in fig. S21B with two local minimums corresponding to vapor and glass, respectively. The local maximum separating the vapor region and glass region (red dashed line in fig. S21B) defines the vapor-liquid interface y 0 . With y s and y 0 identified from fig. S21, A and B, respectively, the liquid-glass interface is y 1 = y 0 + (y 0 − y s ) because the density profile is centrosymmetric.

Different definitions of the surface layers
The three interfaces defined through structural profiles are vapor-dense vapor interface y s , dense vapor-liquid interface y 0 and liquid-glassy layer interface y 1 . In the main text, they are defined as ρ(y = y s,0,1 ) = 5%, 50%, 95%, respectively. Landau theory can provide an alternative definition of

Interface profiles
The density or other structural parameters usually exhibit a centrosymmetric distribution profile in In addition, the surface profile in crystal premelting has been predicted as Eq. S7 from Landau theory (45)

Power law in polymer thin-film glasses
The surface mobile layer is usually studied by comparing the properties of thin-film glasses with different thicknesses (72). For example, for polymer thin-film glasses with low molecular weight, their glass transition temperature T g * and film thickness d can be empirically fitted as follows (72): where T g is the bulk glass transition temperature. The low T * g for a thin film is due to the liquidlike surface mobile layer whose thickness increases with temperature and diverges at bulk T g (72). A thin-film glass can be approximately regarded as the surface region of a bulk glass; thus, the empirical Eq. S12 essentially describes the premelting behavior in Eq. 1 of the main text. However, premelting in thin-film or bulk glasses has not been proposed or compared with the Landau prediction of Eq. S7 for crystal premelting.

Surface melting in ordinary and ultrastable glasses
Most facets of a crystal's surface exhibit surface premelting (15). Such surface liquid can be viewed as a huge postcritical nucleus once the crystal reaches the melting temperature; thus, the crystal melts heterogeneously from free surfaces and preempts liquid nucleation from the interior of the bulk (100).
Consequently, superheated crystals and homogenous crystal melting rarely exist in nature. By contrast, ordinary glasses melt homogeneously within the bulk, and surface melting is negligible because bulk is much larger than the surface region. Ultrastable glass was discovered in 2007 (32) and its surface melting was observed in 2009 (5). The surface melting front propagates to the maximum depth l c where it meets the melted bulk (8,9,16). Therefore, the surface region within l c exhibits surface melting and the deep bulk region exhibits homogeneous bulk melting. This unified picture can explain the observed bulk melting in ordinary glasses and surface melting in ultrastable glasses (8,10): If l c is too small to be observed in a macroscopic experiment, then surface melting cannot be detected. It can be qualitatively expected that the melting behavior depends not only on glass stability but also on heating rate. At sufficiently slow heating, even ordinary glass can melt from the surface, and at sufficiently fast heating, even ultrastable glass can melt from within the bulk.
More quantitatively, the velocity of the melting front v and the bulk glass-to-liquid transformation time τ trans can determine l c as (8): Glass melting experiments (101,102) and simulations (8,12) have revealed the empirical equation where γ is a constant less than 1. The liquid relaxation time τ α follows the Vogel-Fulcher-Tammann (VFT) equation where T 0 is the ideal glass transition temperature. A is a constant. The bulk glass-to-liquid transformation time exhibits the empirical relation (103): Despite their similar forms, Eq. S16 is for the melting transition time for glass, whereas VFT equation (Eq. S15) is for the bulk relaxation time for supercooled liquid.
Eq. S13 is insufficient to identify how l c changes with the heating rate because v increases while τ α decreases with the heating rate. Thus, we further derive the following equation to identify the influence of heating rate on l c on the basis of Eqs. S13-S16 in the literature: A glass increases with the glass stability and A glass > A (103), thus A glass − γA increases with stability, resulting in a larger l c for ultrastable glass than normal glass. Therefore surface melting is more difficult to observe for normal glass. Besides the stability effect, high heating rate causes high transition temperature T (104). Consequently, l c is smaller under high heating rate, and surface melting is more difficult to observe.

Glassy
Bulk      Movie S4. Evolution of log(ρ(y)) and − log(DW (y)) of the monolayer sample under the fast temperature change, which corresponds to Fig. 7D-F.